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Abstract. We start from a study of the density-potential relation for classical 
homeoids in terms of an asymptotic expansion for small deviations from spherical 
symmetry. We then show that such expansion is a useful device that allows us 

5_l ' 

to construct a variety of exact density-potential pairs with spheroidal, toroidal, 
or triaxial shapes for which the deviation from spherical symmetry is finite. As 
concrete analytical applications, we describe: (1) The construction of a family of 
toroidal axisymmetric density-potential pairs one of which is associated with a 
perfectly flat rotation curve (for a member of this family, the supporting two- 
integral phase-space distribution function is obtained in closed form); (2) The 
determination of the aperture velocity dispersion in a wide class of two-integral 
axisymmetric models not stratified on homeoids with central black hole, which 
may be useful for the discussion of the dynamical contributions to the character- 
istics of the Fundamental Plane of early-type galaxies; and (3) For such class of 
models, the construction of the v/cr-ellipticity relation, often considered to assess 
the role of rotation in the structure of elliptical galaxies. 
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1. Introduction 



For the discussion of many astrophysical problems where gravity is important, a major 
difficulty is set by the potential theory. In general, to calculate the gravitational potential 
associated with a given density distribution one has to evaluate a three-dimensional 
integral. Except for special circumstances, where a solution can be found in terms of 
elementary functions, one has to resort to numerical techniques and sophisticated tools, 
such as expansions in orthogonal functions or integral transforms. 

Under spherical symmetry, the density-potential relation can be reduced to a one- 
dimensional integral, while for axisymmetric systems one is left in general with a (usu- 
ally non-trivial) two-dimensional integral. As a result, the majority of explicit density- 
potential pairs refers to spherical symmetry and only a handful of axially symmetric pairs 
are known (e.g., see Binney & Tremaine 1987, hereafter BT87). In few special cases (in 
particular, when a density-potential pair is available in a suitable parametric form) there 
exist systematic procedures to generate new non-trivial density-potential pairs (e.g., see 
the case of Miyamoto & Nagai 1975 and the related Satoh 1980 disks; see also Evans & 
dc Zceuw 1992). 

For non-axisymmetric systems the situation is worse. One class of triaxial density 
distributions for which the potential can be expressed in a tractable integral form is that 
of the stratified homeoids, such as the Ferrers (1887) distributions (e.g., see Pfenniger 
1984, Lanzoni & Ciotti 2003, hereafter LC03) and special cases of the family considered 
by de Zeeuw & Pfenniger (1988). Additional explicit density-potential pairs are given 
by the Evans (1994) models and by those constructed with the Kutuzov-Osipkov (1980) 
method (see also Kutuzov 1998). 

In this paper we draw attention to an elementary yet curious property of the asymp- 
totic expansion for small flattening of the homeoidal potential quadrature formula. Such 
expansion can be traced back to the treatise on geodesy by Sir H. Jeffreys (1970, and 
references therein; see also Hunter 1977). Recently, it has been applied to the modeling 
of gaseous halos in clusters (Lee & Suto 2003, 2004) and to the study of the dynamics 
of elliptical galaxies (Muccione & Ciotti 2003, 2004, hereafter MC03, MC04). The ap- 
parently unnoticed property of the expansion is that it offers a device to construct, in a 
systematic way, density-potential pairs with finite deviations from spherical symmetry. 
In turn, these can be used to carry out a number of calculations explicitly, thus allowing 
for a variety of interesting applications to stellar dynamics. 

The paper is organized as follows. In Sect. 2 we present the basic asymptotic expan- 
sions and we outline the method. In Sect. 3 we study the Jeans equations for axisymmetric 
fluid models associated with density distributions obtained from our method. In Sect. 4 
we focus on three significant applications. In Sect. 5 we make some concluding remarks. 

Send offprint requests to: L. Ciotti 



L. Ciotti and G. Bertin: Exact density-potential pairs 3 

In Appendix A we address the issue of the relevant parameter space. In Appendix B 
we clarify the relation between the homeoidal expansion and the expansion in spherical 
harmonics. In Appendix C we summarize the intrinsic kinematical profiles of a family of 
axisymmetric galaxy models with central black hole. 




Fig. 1. Isodensity (top), constant projected density (middle), and isopotential (bottom) 
contours (in arbitrary units) for the toroidal model described by Eq. IjlSJ) with a = 3.1 
(left) and a — 4.9 (right). The lower value of a is near the critical value for which 
the surface brightness is everywhere infinite, while the upper value is near the critical 
value for which the central mass diverges. The coordinates R, y, and z are normalized to 
the scale-length a. Note how a strongly non-spherical density distribution, such as that 
represented on the top right, produces a nearly spherical potential (bottom right). 
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2. The method 

We start by assuming a density distribution stratified on homeoidal surfaces, p(x) = 
pop(m), where x = (x, y, z), po is a density scale, and 

2 ? 2 2 2 9 

2 x ^_ y ^ z x ^ y • ^ z 

a 2 b 2 c 2 a 2 a 2 (l — e) 2 a 2 (l — T]) 2 ' 

with a > b > c > 0, b/a = 1 — e, and c/a = 1 — t]. The prolate case corresponds to 
e = 77 < 1 and the oblate case to e = and < 77 < 1. In models of astrophysical interest, 
the density p may extend to infinity (with finite or infinite mass) or be truncated. For 
simplicity, in this paper we discuss only non-truncated density distributions. 

It is well known (e.g., see Kellogg 1953, Chandrasekhar 1969) that the (inner) poten- 
tial associated with a homeoidal density distribution can be written as 
r°° A$(x- e ??) 

<Kx) = -irabcpoG / ^ ' ' " dr, (2) 

Jo A ( T ) 



where A(t) = y/ (a 2 + t) (b 2 + t) (c 2 + r) , = 2 f™r x - T ) p(m)m dm, and m 2 (x;r) = 
x 2 1 (a 2 +T)+y 2 1 (b 2 +r)+z 2 / (c 2 +t). If we rescale density and potential to the quantities po 
and 4irGpoa 2 , respectively, the Poisson equation for the dimensionless density-potential 
pair (p, <p) becomes V 2 = p where V 2 = d 2 /dx 2 + d 2 /dy 2 + d 2 /dz 2 , i = x/a, y = y/a, 
z = z/a and 

^=-(l-e)(l-v)-J Q A(r) dr. (3) 
This will be the starting point of our analysis. We now expand the density p(m) and the 
potential 4> given in Eq. J2J up to the first significant order in the flattening parameters 
and find: 



,>[ m ) = f,{ r) + ey2 + ^ 2 p\f) + Q(e 2 + if), 1 -1 1 



r 

where p' = dp(m) / dm evaluated at e = rj = 0, and 

4> = Mr) + + \Mr) - <h(f)] + ( £ y 2 + ^ 2 )0 2 (r) + 0(e 2 + rj 2 ), (5) 



where 



Mr) 



-— it I p{m)m A dm — — / p(m)m dm, (i = 1); (6) 



p{m)m dm— j p(m)m dm, (i = 0); 
f 

1 r 

— I p{m)m dm, (i = 2). 
r ° Jo 

In the above expressions, the dimensionless spherical radius is defined as f = r/a 



y/x 2 + y 2 + z 2 . Note that ^o(^) is the potential associated with /5(f), as can be seen 
from its explicit expression. The infinite upper limit of integration appearing in two 
integrals could be changed to a finite constant if desired; this property is useful to deal 
with density distributions for which the two integrals would diverge. Similar formulae 
can be obtained for the expansion at fixed total mass of density distributions of finite 
total mass (e.g., see MC04, Eq. [C.l]). 
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2.1. Generating new families of non-spherical density-potential pairs 

By ordering arguments and by the linearity of the Poisson equation, it follows that 
the truncation of Eqs. iQ-jSj to first order in the flattening parameters produces exact 
density-potential pairs independently of the value of e and 77. In order to be physically 
acceptable, the density distribution truncated to first order in e and 77 must be posi- 
tive definite (see Appendix A). Density-potential pairs belonging to the families thus 
constructed can be non-spherical at the "non-linear" level, because the values of the flat- 
tening parameters can differ significantly from zero. On the other hand, the truncated 
potential is also the approximate potential for the density p — pop~(m) when e — > and 
77 -> 0. 

Now we proceed one step further and use Eqs. (0}-© to construct new density- 
potential pairs, independent of the ellipticities e and rj defining the homeoidal expansion. 
In fact, by the same ordering argument used above, it is immediate to see that </>i(r) — 
4>o(r) + z 2 (f)2(r) is the potential associated with the "density" z 2 p'(r)/f. For radially 
declining density distributions /3(f), we can thus consider the (dimensionless) pair 

g = ^lM ; § = ^-^-^. (7) 
r 

Two similar density-potential pairs can be then obtained trivially by renaming variables, 
i.e., by replacing z with y or x in Eq. JjJ. We now combine the three density-potential 
pairs linearly to get the following generally "triaxial" solution 

£=(ai 2 +/3y 2 + 7 5 2 )^P, (8) 

r 

$ = (a + 13 + 7) 0o (f) - 0i (f )] - (ax 2 + I3y 2 + 7 z 2 )<fe (f) . (9) 

If we now set a = (3 = 1 and 7 = 0, we find a new axisymmetric family characterized by 
toroidal geometry 

Q = i?\P^A. $ = 2(0 o -0i)-i? 2 02. (10) 



Here R = \/x 2 + y 2 is the dimensionless cylindrical radius. 

That Eqs. (|7|l- l|lU|l indeed identify explicit density-potential pairs can be checked 
a posteriori by deriving the potential $ from the density g by means of a standard 
expansion in spherical harmonics (Appendix B). Note however that, if we started with a 
seemingly simpler idea, i.e., to transform the density p(m) in spherical coordinates and 
then to expand it in spherical harmonics, we would have discarded all the expansion 
terms above the monopole as candidates for acceptable density distributions because the 
related mass vanishes. In contrast, the present method works because, at least to first 
order, it is possible to isolate a density term with constant (in our case, negative) sign 
over all the space. 

In the approach followed so far, we have considered the density p — pop~(m) as given: in 
this way, several density-potential pairs can be easily obtained, for example by starting 
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from well known seeds for /5(f), such as the density of the 7-models (Dehnen 1993, 
Tremaine et al. 1994, see also Appendix A). However, the generated density profiles 
given in Eqs. (|7|)- (|10fl remain largely out of control. 

On the other hand, it is possible to specify the spherically symmetric factor p = 
\p'(f)\/r in Eqs. 0, JHl, and (|10f) : the potentials 4>o, 4>i, and cj>2 required to complete the 
non-spherical density-potential pair are recovered by inserting in Eq. JSJ the supporting 
density 

p{f) = J_ P(t)tdt. (11) 

The function p should go to zero sufficiently fast when r — > 00, in order for the integrals 
appearing in Eq. (fTTf) to converge. We will use this approach to build the special model 
presented in Sect. 4.1. 

Another way to specify the supporting density p (and so p) is by imposing a given 
rotation curve in the equatorial plane, that is a circular velocity v c {R) defined from the 
relation v 2 (R) = v 2 (R) / (AnG poa 2 ) — Rd^/dR at z = 0. In fact, by considering the 
dimcnsionless potential $ in Eq. (10) and by assigning the function v 2 , one gets and can 
solve the following (inhomogencous Euler) differential equation for p: 

1 d 



2 ^ 

3-2 



dR 2 dR dR 



R dR 



0. (12) 



3. The Jeans equations of stellar dynamics for axisymmetric 

density-potential pairs (embedded, if desired, in a spherical external field) 

In this Section, in view of applications to galactic dynamics, we study the Jeans equations 
for axisymmetric density-potential pairs either in the form of Eqs. J7J and i|10|l . or in the 
axisymmetric cases of Eqs. In our study we can also include the case in which 

the density distributions mentioned above are added to a separate spherically symmetric 
density distribution. 

In all the axisymmetric cases described above, the pair (g, 3>) belongs to the general 
family of density-potential pairs 1 

P[R] = A(r) + R 2 B(r), cj> [R] = C(r) + R 2 D(r), (13) 

where the subscript [R] means that the pair is made explicit with respect to R by using 
(if required) the identity z 2 = r 2 — R 2 . The possibility to switch to an alternative repre- 
sentation (i.e., from p^ and ^rra to pr 2 i and 4>\ z ]) is useful for the solution of the Jeans 
equations (when no subscript appears, it means that either form can be considered). 
Note that for D{r) cx 1/r 4 the potential belongs to the class of Stackel potential (Stackel 
1890). 

Let us assume that the density-potential pair defines the equilibrium configuration of 
a stellar dynamical system for which the underlying distribution function (DF) depends 
1 In the self-gravitating case A(r) oc C"(r) + 2C'(r)/r + 4D(r) and B(r) oc D"(r) + 6D'(r)/r. 
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7 



on the two classical integrals E and J z (star energy and z-component of the angular 
momentum, per unit mass), so that the velocity dispersion tensor <7y is diagonal and 
a 2 = <jrr = cr zz ; we will denote by crH, the dispersion a w . Along the azimuthal direction, 
part of the kinetic energy can be stored in systematic motions 2 , so that v 2 = u 2 + er 2 . 

Given the form of the density-potential pair in Eq. (|13|) . the Jeans equations (e.g., 
see LC03, Eqs. [7] -[8]) can be easily solved as 



Clearly an ambiguity remains as to the amount of support related to systematic motions 
with respect to that of random motions in the toroidal direction. A commonly used 
decomposition assumes m 2 = fc 2 (v 2 — a 2 ), and so a 2 = a 2 + (1 — fc 2 )(v 2 — a 2 ), where k 
is generally taken to be a given constant (Satoh 1980, but see Ciotti & Pellegrini 1996). 
When k = the system is azimuthally supported only by velocity dispersion, and no net 
rotation is present; in the opposite limit, k = 1, the velocity dispersion is isotropic. Of 
course, the above procedure can be applied only when v 2 — a 2 > everywhere. 

The general expression of projection integrals for axisymmetric systems can be found 
in LC03. For simplicity, in this paper we focus on the case when the system is observed 
edge-on (i.e. from a line-of-sight where the contribution of rotation is expected to be 
largest), even though for density distributions as in Eq. 1|13[) explicit projection formulae 
can be obtained for any inclination of the line-of-sight direction n. Thus, we project 
along the x axis, so that n = (1,0,0) and the projection plane is the (y,z) plane. The 
projected density is given by E = pdx, while the mass-weighted projected streaming 
velocity field is given by Eu p = pu n dx, where u n = u^n,; = u x = —u^sinip and 
sin ip — y/R; similarly, the projection of the square of the streaming velocity field is 
T,V 2 = pu^dx. Finally, the line-of-sight velocity dispersion field is af os = a 2 +V 2 —u 2 , 
where Ecr 2 = pcr 2 dx, and a 2 = aijUiUj — a xx = a 2 + (1 — fc 2 )(v 2 — a 2 ) sin 2 ip (see 
Eqs. [B5] in LC03). 

3.1. The problem of the supporting distribution function 

At the end of this Section, we may ask whether it would be possible to reconstruct the 
phase-space DF for the general family of densities given in Eq. I|10[l . As is well known, 




(14) 



and 




(15) 



the DF 



/ = f J 2 z h(£)O(£), 



(16) 
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where £ = —v 2 /2 + <3> is the star binding energy per unit mass, and is the Heaviside 
step function, produces a mass density given by 



without streaming velocity (e.g., Hunter 1977; for a more general family of distribution 
functions containing Eq. i|16|l as a special case, see also Ciotti, Bertin & Londrillo 2004). 
Equation (17) can be used to determine the two-integral DF associated with a density 
profile of the form given by Eq. I|10|) . in the limiting case when the total potential — 
is dominated by an external, spherically symmetric potential — "^cxtM (so that, in such 
approximation, we ignore the contribution to the potential given by the field produced by 
g itself). In particular, three different strategies can be considered. In the first two cases 
h(£) is assigned and then Eq. (17) is used to determine either the spherical factor p(r) for 
given v I / ext(^), or, alternatively, \& ex t(r) for given p(r). In the third case, the two functions 
^ext(?") and p(r) are assigned, and h(£) is obtained by standard Abel inversion of the 
function p[ v I' 0X t]. Of course, in each of these cases the appropriate physical requirements 
must be satisfied, namely the mass density distribution associated with — ^cxt via the 
Poisson equation and the function h must be non-negative. Note that if dp/dr < 
then r provides the natural parameter to obtain the mapping <f> = Q(g, R 2 ) from Eq. 
I)10[l. This is an important step, because g — g(Q,R 2 ) is the function needed to recover 
the two-integral DF of axisymmetric systems (cf. the Fricke [1952] series expansion, the 
Lynden-Bell [1962], the Dejonghe [1986], and the Hunter [1975] integral transforms, and 
finally the Hunter & Qian [1993] integral contour method; see also Sect. 4.1.2). 

4. Three applications 

Following the normalization introduced at the beginning of Sect. 2, in this Section and 
in Appendix C, the density g, the relevant radii, and the potential $ are meant to be 
normalized to po, to the scale- length a, and to AwGpoa 2 , respectively. It follows that 
masses, projected densities, and velocities are normalized to poa 3 , poa, and 2a^/nGpo; 
for simplicity, we omit the tilde above these quantities. 

4.1. Power-law tori 

In this subsection we describe a toroidal stellar system belonging to the family associated 
with Eq. I|10l) . for which the general methods outlined in Sect. 3 allow us to obtain 
explicitly the relevant kinematical profiles (and their projections on the plane of the 
sky), even in the presence of a central black hole. For simplicity, we focus here on the 
fully self-gravitating case, i.e. on the case when no black hole is present. The simplest 
model is that of a power-law torus, defined by the (normalized) density distribution 




(17) 



Q = 



— , (a>0). 



(18) 
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The total mass of the model is infinite, independently of the value of a, an unpleas- 
ant property shared with spherically symmetric (non-truncated) power-law models. The 
central mass is finite for a < 5. According to Eq. the supporting density is 

p = r 2 ~ a /(a - 2) (for a > 2), and so from Eqs. {TUJ) and © 

4 f2 



$ = 



(a-2)(7-a) 



- ( 21nf - -^r 




(19) 



For 4 < a < 5 the gravitational potential vanishes for r —* oo and diverges at the 
origin: thus, only orbits with negative total energy are bound. For 2 < a < 4 an upper 
truncation must be applied to some integrals appearing in Eq. 10; the additive constant 
resulting in Eq. (|19|l has been set equal to zero, and the potential vanishes at the origin. 
Finally, for a = 4 the potential diverges both at r — and r — oo and so orbits are 
bound independently of their energy. 

For 2 < a < 5 and a^O the circular velocity in the equatorial plane is 
2 _ {-a 2 +9a-l6)R A - a 



c (a-2)(5-a)(7-a)' 
while v 2 = 2/3 for a = 4. Curiously, the circular velocity vanishes for a — 



(20) 



9/2 — VT7/2 ~ 2.44. In fact, for 2 < a < a n Eq. (|20|l would predict a negative square 
circular velocity, in the sense that in this range of a the radial force field in the equatorial 
plane, is directed outward. For 2 < a < 5, Eqs. 114|l - l|15Jl give 



go - 
and 
g(v 2 



ft2 p2(l-a) 

7 — a 



2f 2 



R 2 



(a- 2) 2 (5 -a) 2(a - 1) 



2 ^2~2(l-q) 

(a-2)(7-a) 



2r 



,-.2 



I? 2 



(a-2)(5 



(21) 



(22) 



since the r.h.s. of Eq. I|22l) is positive, the Satoh decomposition can be applied. 

In the isotropic case (k = 1), Eq. (1221 gives the square of the streaming velocity it 2 ; 
in this case in the equatorial plane (where r = R), the ratios a 2 /v 2 , v^/v 2 , and u^/v 2 are 
constant. The isotropic case can also be interpreted as the description of a toroidal fluid 
structure for which the velocity u v is not constant over cylinders (corresponding to the 
fact that the system is baroclinic; e.g., see Tassoul 1978): the asymmetric drift v 2 — u 2 v 
on the equatorial plane presented by all the models scales as R 4 ~ a , vanishes for a = 3, 
and is negative for lower values of a. Along the z axis (where R = 0), according to Eq. 
I|21l) the pressure vanishes, and so does the quantity g{v^ — a 2 ): while this latter behavior 
is a common property of two-integral axisymmetric systems (e.g., see Ciotti & Pellegrini 
1996, LC03), the former is a peculiarity of the present toroidal model. Along the major 
axis in the edge-on projection plane (i.e., along the y axis, at z = 0), u p oc ky 2 ~ a l 2 , while 



V 2 



where the proportionality constants are simple functions of 
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a that can be easily computed for a > 7/2; below this value of a the projected velocity 
dispersion diverges. The properties of two toroidal models associated with Eq. I|18(l are 
illustrated in Fig. 1. 

4.1.1. The torus with flat rotation curve and the "neutral" torus 

From the solution of Eq. I|12|l with v c = constant, one can prove that the a = 4 power- 
law torus is the only model of this class characterized by a perfectly flat rotation curve. 
For a generic value of k, along the major axis in the edge-on projection, we have u 2 = 
8fc 2 /(9^ 2 ), and af os = 1/3 - 8fc 2 /(9^ 2 ). 

A more general family of power-law tori with flat rotation curves can be obtained 
by expansion of the oblate spheroid p cx 1/m 2 (the so-called "isothermal" spheroid) 
to higher orders in the flattening, beyond Eqs. In fact, such oblate spheroid is 

characterized by a flat rotation curve on the equatorial plane (for any flattening; e.g., see 
BT87, Chap. 2). Therefore, following the arguments presented in Sect. 2.1 we find that 
the tori 

i? 2n 

£ = ^T2> O integer), (23) 

all belong to the list of density distributions with perfectly flat rotation curve, which 
includes the singular isothermal sphere and the Mestel (1963) disk (see also Monet, 
Richstone & Schechtcr 1981; the toroidal solutions presented by Toomre 1982 can also 
be expressed as a superposition of power-law tori). When projected face-on, much like 
the isothermal sphere, the family given in Eq. 11' Ml) produces the l/R surface density 
distribution of the Mestel disk. 

A torus with very unusual characteristics is the following. If we discard one of the 
two homogeneous solutions of Eq. Ijl2|l (corresponding to v c = 0), because of its infinite 
central mass, we are left with the solution p cx f- 5 / 2 +vT7/ 2 . From Eqs. (20) and (|11(1 
we see that the related density is given as in Eq. |JTSj, with a = a n . This torus is thus 
characterized by a constant potential in the equatorial plane, so that the circular velocity 
vanishes, as already pointed out in the previous Section. From a mathematical point of 
view, adding this neutral torus to a disk, would not alter the rotation curve in the disk. 
In other words, if we add the neutral torus to a flat-rotation density distribution, the 
resulting density distribution will be associated with a flat rotation curve of the same 
amplitude as the original. 

We conclude by noting that fully self-consistent, self-gravitating toroidal stellar sys- 
tems corresponding to the density-potential pair (|18fl - i|19[) with a < a n do not exist. In 
fact, for 2 < a < a n the zero- velocity curves in the meridional plane (i?, z) are open. 
Surprisingly, the Jeans equations would lead to apparently "innocent" solutions even 
for a < a n . This concrete example illustrates the risks that are taken when stellar sys- 
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terns are modeled by means of the Jeans equations without a final check on the physical 
viability of the solutions found from the formal analysis. 



4.1.2. The two-integral distribution function of power-law tori 

The general considerations presented in Sect. 3.1 can be elaborated in greater detail 
for the family of power-law tori. In the case of non-self gravitating power-law tori, the 
distribution function in the presence of a dominant massive central black hole of mass 
Mbh can be explicitly written as a Fricke distribution function by adopting h{£ ) = 
£«-5/2 ( with a > 3/2) an d f Q = 3 / 9 a a - 2 /[2 7 / 2 7r(GM B H) Q B(5/2, a - 3/2)] in Eq.JHJ, 
where B(x,?/) is the standard Euler complete Beta function. 

As noted at the end of Sect. 3.1, by combining Eqs. (|18fl and i|19[) . the function 
$ = &(g,R 2 ) can be obtained explicitly. In addition, for the entire family of power-law 
tori (with a > a n ), the envelope and the dynamical window (relevant to such method) 
can be expressed in terms of elementary functions: for example, the envelope is defined 
by the relation J 2 oc exp(3i? c ) when a = 4, while J 2 cx lE^I 3 ^ in all the other cases, 
where E c is the energy per unit mass of a star orbiting on the circular orbit. Finally, for 
two members of this family it is possible to construct the function g — g{^, R 2 ) explicitly, 
thus providing analytically all the required ingredients for a numerical implementation 
of the Hunter-Qian method (Qian et al. 1995). In particular, 



(^4$2 + 2^2 _ 2 $)S 
R 1 

4W 2 (i? 2 e- 3< V2) 
R 2 



(a = 3), 



(24) 



(a = 4), 



where W(x) is the Lambert function 3 . Remarkably, for a = 3 the density g can be 
expanded in series of ip 2 = i? 2 /(2<I> 2 ), and the associated series of Fricke's terms can be 
re-summed in closed form as 



f{E,J z ) 



945 £ 

211/2^7^72" 

k sign( J z 



128tt E 5 / 2 



G 

4F3 



6 ' 4y 256 
3 5 11 13 5 7 



143? ^ (17 19 r 27? 
2^1 r- 



6 ' 6 



4 _27? 
2' 2' 6 ' 6 ' 4' 4' ' 4 



(25) 



which is positive for < ? = J 2 /(4E 3 ) < 32/27 and < k < 1 (where the upper 
limit refers to circular orbits and p F q is the standard hypergeometric function, e.g. see 
Gradshteyn & Ryzhik 1980). Note that the odd part has been recovered under the Satoh 
ansatz, for which u v — kyf<$> from Eqs. (18), (19), and (22). We verified by direct inte- 



gration over velocity space that the DF returns indeed the a — 3 torus. 



3 The Lambert function is defined by the identity lnW(e z ) + W(e z ) = z. 
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4.2. Scale-free oblate axisym metric systems treated in terms of the Jeans equations 
and the Satoh decomposition 

We now apply the formulae derived in Sect. 2.1 to the family of power-law oblate (e = 0) 
spheroids. Thus, we start from p = py/m 7 , where the range < 7 < 3 is such that the 
density profile is radially decreasing and the central mass remains finite. 

Several properties of these models have been investigated by Qian et al. (1995), who 
reconstructed their phase-space properties; in addition, Evans (1994) and Evans & de 
Zeeuw (1994) studied the observational properties that would result from a decomposition 
similar (but not identical) to the Satoh (1980) decomposition for a similar class of models 
(which they called "power-law" models, in their scale-free limit). 

Here we study the relevant Jeans equations on the basis of the Satoh decomposition, 
which has found widespread applications to elliptical galaxies (e.g., see Binney, Davies 
& Illingworth 1990; van der Marel, Binney, & Davies 1990; van der Marel 1991). We are 
then able to provide expressions for the intrinsic and projected kinematical profiles in 
a form that may be useful to model the observations. These expressions can be easily 
generalized to the case in which a central black hole or a power-law dark matter halo 
with different flattening are present. 

From Eqs. ©-© we have 



where the arguments described in Appendix A impose 77 < I/7; the limit r\ = 1/7 recovers 
the case described in Eqs. CHI- El- The presence of a central black hole of mass Mbh 
can be modeled in Eq. (|27|l by adding a term — pbh/^, where /ibh = MBn/^Po a3 - The 
density in Eq. (|26l) becomes prominently peanut-shaped at high values of 77; this is a 
common property of the original model studied by Schwarzschild (1979), of the flattened 
Plummer sphere (Lynden-Bell 1962), and of the family of "power- law" models (Evans 
1994; the family includes the Binney [1991] logarithmic potential as a special case). 

For 7 > 2 the potential vanishes at infinity and is negative divergent at the origin. For 
< 7 < 2 an upper truncation must be applied to some integrals appearing in Eq. 
in Eq. 1271) an additive constant has then been set to zero, so that the potential diverges 
at infinity and vanishes at the origin. Finally, for 7 = 2 the model is characterized by a 
flat rotation curve, as the seed density distribution. In fact, 



1 — 777 jrjR 2 

fi pr+2 ' 



(26) 




(27) 




5 — 7 — 2?7 



+ 



Mbh 
R 



(28) 



(5-7)(3-7)^" 2 
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In the entire range < 7 < 3, the models do not exhibit the repulsive behavior noted at 
the end of Sect. 4.1.1 for power-law tori and so, from this point of view, are physically 
acceptable. 

The fully nonlinear (i.e., up to rf terms included) Jeans equations (f 1 4|> - (| 1 5|> for the 
assumed density-potential pair can be solved without difficulty. Here, for simplicity, 
we use the solution up to first order in rj; the analytical expressions of the intrinsic 
and (edge-on) projected kinematical fields are given in Appendix C. In Appendix C we 
also show that the projection of all the quadratic (intrinsic) kinematical profiles can be 
obtained explicitly even for /kbh 7^ 0. Unfortunately, the streaming velocity field cannot 
be obtained in explicit form when /xbh 7^ 0; for this case we can only provide asymptotic 
expressions for small values of the projected radius. Close to the origin, independently 
of the value of 7, the dynamical effects of ^bh are always dominant. 

For /iBH = 0, the analysis presented in Appendix C gives (for 3/2 < 7 < 5/2) 



M P a 2 a = 



7r 3 / 2 r(7-3/2)i 5 - 2 ^ 

2(5- 7 )(3-7)(5-2 7 )r( 7 ) 



2(5 - 7) + r]{2Y - 7 7 - 2) 



2 „ n , 4fcVW4-i/2)r( 7 /2)r( 7 ) 



r2(3 7 /4)r( 7 /2-l/2)r( 7 -3/2) 



(29) 



where I = yjy 2 + z 2 is the (dimensionless) circular aperture radius and T is the standard 
Euler gamma function. 

For /iBH 7^ and I — > (that is, in practice, for apertures not larger than the sphere 
of influence of the black hole) 

„ 2 /W/ 2 r( 7 /2 + i)^ ^ 



2 7 (2 - 7 )r( 7 /2 + 5/2) 
2(3 + 7) -7(4-7)77 - 



fc 2 r ?7 (7 2 -l)r 2 (7 + l/4) 



(30) 



r 2 ( 7 /2 + 3/4) 

for 1 < 7 < 2. Thus, scale-free oblate systems with a central black hole do have finite 
projected fields in the range 3/2 < 7 < 2. In Eqs. (29)-(30) c 2 is the aperture velocity 
dispersion defined as <r 2 (^) = J y 2 +Z 2 <t >2 ^o~i oa dydz / 'M p (£) , where M v {£) is the projected 
mass within I. 

By combining the two formulae above we obtain an expression for the aperture ve- 
locity dispersion which is asymptotically correct for small and large apertures when 
Mesh 7^ 0. 



4.2.1. Effects of rotation on the aperture velocity dispersion 

In LC03, the problem of projection effects on the tilt and thickness of the Fundamental 
Plane (FP) of elliptical galaxies was investigated by means of a novel approach, intro- 
duced by Bertin, Ciotti & Del Principe (2002). In particular, it was shown that rotation 
has no significant effects on the value of the measured "central" velocity dispersion (con- 
sistent with the results obtained by van Albada, Bertin & Stiavelli (1995) for cuspier 
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models), and thus has a negligible effect on the FP tilt and thickness. However, the 
models used suffered from possessing unrealistically flat density profiles in their central 
regions. Here we illustrate how such legitimate concern can be resolved with our method, 
for simplicity restricting our discussion to the edge-on case, where the effect of rotation 
is maximum. 




1 .24 lo0(T o I 0.82 log T e I 0. 1 82 

Fig. 2. Relative displacement of galaxy models induced by different amounts of rotational 
support with respect to the edge-on FP for different slopes 7 of the density profile and 
different flattenings. The two dotted lines parallel to the FP line represent the 1-rms 
thickness of the FP. The aperture velocity dispersion of each model is normalized to 
that of the corresponding isotropic rotator (k = 1), which is placed on the FP (solid 
diagonal line). Dotted and short-dashed lines (labelled by 7, lower part of the figure) 
represent models of Sect. 4.2 without and with a central dominant black hole, in the case 
of maximum flattening, except for the dotted line labeled (2.4, 1/4.8) which represents a 
model with smaller flattening. Long-dashed and dotted-dashed lines refer to Eq. (Cll) 
of LC03, i.e., to the velocity dispersion of the homogeneous oblate ellipsoid (n = 0) 
calculated within an elliptical aperture with semimajor axis equal to the effective radius 
R e (label "e") and over the whole model (label "t"), for two different flattenings (axial 
ratio q = 0.3 and q — 0.6). 
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In Fig. |2 we show the position of our models with respect to the edge-on FP. The 
models are first (arbitrarily) placed on the FP (the diagonal solid line) in the isotropic 
case (i.e., k = 1); the two dotted lines parallel to the diagonal show the 1-rms dispersion 
of real galaxies around the FP. We refer to the FP coefficients measured by J0rgensen, 
Franx & Kjaergaard (1996; for Coma cluster ellipticals and Hq = 50 km s _1 Mpc^ 1 ). 

The models are then "moved" horizontally by letting their "observed" velocity dis- 
persion ct vary with the Satoh parameter k down to k = 0, while their scale-length and 
mean surface brightness are kept constant. To better illustrate our point, in Fig. [21 we 
place different models at different values of R c . 

Because of the scale-free nature of the models considered the relative effects are 
independent of the adopted aperture (cfr. Eqs. [29]- [30]). In all models oo increases for 
decreasing k, i.e., as the flattening becomes more and more supported by the azimuthal 
velocity dispersion. However, the displacement depends on the density slope 7 and on 
the model flattening. Note how cuspier models are more affected by rotation than softer 
ones. The effect of flattening at fixed 7 is then shown by the line labeled (2.4, 1/4.8): for 
this model the density slope is 7 = 2.4, but the flattening is only half of the maximum 
allowed. The effect of rotation for models with a dominant central black hole is apparently 
weaker. 

The behavior of cuspy models, even in the presence of a massive central black hole, is 
not very different from that of the Ferrers models in LC03, thus confirming the findings 
of van Albada et al. (1995) and the statistical results of LC03. 

4.3. The v/a-ellipticity plane 

One important property of a family of models characterized by different degrees of flat- 
tening, is the way it covers the so-called v/cr-ellipticity plane. In fact, such plane is a 
simple tool often used to investigate to what extent the observed flattening of galaxies 
should be ascribed to rotation (e.g., Illingworth 1977, Binney 1978). For the purpose, the 
observed quantities [v/cr, a measurement of the ratio of observed rotation to velocity dis- 
persion, and the observed ellipticity) are usually compared to a curve associated with the 
behavior of classical spheroids (see Chandrasekhar 1969). However, such a comparison 
may be misleading, because (1) galaxies need not conform to the properties of classical 
spheroids with density stratified on homeoidal surfaces (Roberts 1962), and (2) the ob- 
servations sample the central regions and may thus provide insufficient information for a 
proper comparison with the theoretical expectations. Indeed, Evans & de Zeeuw (1994), 
in the study of their "power-law" models, noted that their nearly isotropic models are 
associated with points in the v/d-ellipticity plane systematically below the curve of the 
classical spheroids, which is obtained from virial (and therefore integrated) quantities. 
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Fig. 3. The v/cr-ellipticity plane for edge-on isotropic rotators (BT87, Eq. [4.95]). Here 
the ellipticity is defined as the value of tj. The thick solid curve is the locus of classical 
(oblate) spheroids. The dotted line "c" refers to the n — model in LC03, where the 
projected rotational velocity and the aperture velocity dispersion are taken at the "core" 
radius and within an elliptical aperture of semimajor axis equal to the core radius, re- 
spectively. The dotted line "0" refers to the same model, where the projected rotational 
velocity is taken at ^/2/5i? t (see text), while the aperture velocity dispersion is measured 
at the center. Long-dashed and short-dashed lines refer to our models, with and without 
a dominant central black hole respectively, and for different values of the density slope 
7 (given by the labels) of po/m 7 . 



Here we address this issue by studying the behavior of the family of models introduced 
in Sect. 4.2. Note that, even though our models are indeed obtained from a seed density 
distribution stratified on homeoidal surfaces, their density distribution is not stratified 
on homeoids, being the sum of a spherical density distribution and a (negative) toroidal 
density. To stay closer to the issue of the relevant observed quantities we define a in 
terms of the aperture velocity dispersion (Eqs. [29]- [30]) and v in terms of the major 
axis projected streaming velocity (obtained from Eq. [C8]). Note that, because of the 
scale-free nature of the present models, the rotational parameter is independent of the 
aperture radius adopted, when both the streaming velocity and the aperture velocity 
dispersion are measured at the same radius. 
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In Fig. |21 we show the model position in the ellipticity- rotational parameter (v /a) 
space for edge-on, isotropic stellar systems. For the n = model described in LC03 the 
major axis projected rotational velocity increases linearly with radius and thus, following 
Evans & de Zeeuw (1994) (footnote 2 there), we adopt as reference value of v the value u p 
measured at the "core" radius (i.e. at ^/3R t /2, where Rt is the ellipsoid semimajor axis). 
The aperture velocity dispersion a referred to an elliptical aperture of core radius (curve 
"c"). It is thus clear that the rotational parameter vja can stay above or below the locus 
of classical spheroids. In fact, it can be easily proved that for the constant density spheroid 
the line-of-sight velocity dispersion and the projected velocity dispersion coincide (Eq. 
[C.10] in LC03), and so from the projected virial theorem (e.g. Ciotti 1994) the aperture 
velocity dispersion over the whole object coincides with the virial velocity dispersion (in 
the isotropic case). Also, the streaming velocity, projected for the edge-on isotropic case, 
evaluated along the major axis at \j2jhRt coincides with the rotational velocity obtained 
from the virial ordered kinetic energy. For example, curve "0" refers to the n = model 
where v is taken at W2/5Rt, while a is the central projected velocity dispersion. 

The result we wish to point out here is associated with the dashed lines, which 
correspond to our models with (long-dashed) and without (short-dashed) a dominant 
central black hole, for density slopes near the limits of the acceptable range (labels near 
the curves give the value of 7 for the assumed p = po/ra 1 of Sect. 4.2). 

Thus, while our analysis confirms the importance of aperture effects even for the 
classical spheroids, we are also able to show that in general models different from classical 
spheroids may be found well below the classical locus, with deviations increasing at larger 
ellipticities. 



5. Concluding remarks 

The applications addressed in some detail in Sect. 4 are only a small sample of the 
many applications that can grow out of the general method introduced in this paper. 
Among other applications not addressed here, we may mention the following. MC04 
adopt the method to calculate stellar orbits in triaxial elliptical galaxies, while Barnabe 
et al. (2005) investigate the properties of exploratory models for the low-rotation, extra- 
planar gas in spiral galaxies. Finally, we should recall that the interest in the study of 
toroidal structures is also motivated by the fact that centrally depressed surface bright- 
ness profiles need not always be associated with the projection of spherically symmetric, 
non-monotonic density profiles (Ciotti 2000), but may trace instead the presence of truly 
toroidal structures. 
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Appendix A: Limits on e and 77 

When considering homeoidal expansions we have to impose that the density distribution 
thus obtained be positive-definite, which then sets limits on the maximum values of e 
and rj as a function of the seed density profile. Without loss of generality we take e < rj. 
In the case of the expansion at fixed total mass (see MC04, Eq. [C.l]) we require 

ey 2 + r\z 2 d\np(r) 



1 



din'/ 



> 0. 



(A.l) 



If we assume that p(r) is a non-increasing function of r, by defining Am 
sup[ |dln p(r) / 'd\nr\ > 0, the inequality above is then satisfied when 



e> (A M - l)r? - 1, 



(A.2) 



because the factor (ey 2 + r\z 2 )jr 2 attains its maximum value rj on the z axis. Clearly, 
if the seed density profiles are too steep little room is left for the expansion procedure. 
In the unconstrained expansion (Eq. 0]), a similar analysis for e < r\ < 1 leads to the 
condition 1 > AmV- 

To illustrate these conclusions, let us consider the density profile 

' (A.3) 



m a (l + m b ) c ' 

with a > 0, b > 0, c > 0, for which Am = a + be. Thus, in the unconstrained expansion 
we have to require 77 < 1/ (a + bc), while in the oblate case at fixed total mass the relevant 
condition is 77 < l/(a + be — 1). Thus, for a power-law model (c = 0), < 77 < 1/a. 

A particularly interesting case of expansion at fixed total mass is represented by the 
so-called 7- models (Dehnen 1993, Tremaine et al. 1994), for which a = 7, b = 1, c = 4— 7. 
The 7-models received much attention not only for their analytical simplicity, but also 
because the projected density distribution of the 7 = 2 model (Jaffe 1983) is well fitted 
by the de Vaucouleurs (1948) i? 1 / 4 law. Curiously, the positivity condition for the oblate 
case is < 77 < 1/3, i.e, the expansion must be limited to models rounder than E7. For 
the 7 models we have M = 1/(3 — 7), and so, from Eq. (C.l) of MC04 

1 + e + f] (ey 2 + rjz 2 )(4:f + 7) 



Q 



r~<(l + ?) A -~< f 2 +T(l + f) 5 -T 



(A.4) 



with 



1 



(2- 7 )(3- 7 ) 
r 



1 



2-7 



- 1 



(7^2); 



In 



1 



(7 = 2) 
Tf 2 -T(3-7 + f) 



3(2- 7 )(3- 7 ) 
r-2 



3r 2 



1 



(l + f)3- 7 

21n(l + f) 

3f3 : 



f 2 -TF 21 (4- 7 ,5-7,6-7,- 



3(5-7) 



(A.5) 



(7^2); 



(7 



(A.6) 
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i^l(4-7,5-7,6-7,-f) 

' ^ ^ (A 7) 

2 + f 21n(l + f) 1 '> 

f*(l + f) f5 » 17- )• 

The case 7 = 1 was adopted by MC03 and MC04, and similar formulae can also be 
obtained for the so-called modified Jaffe profile, p oc l/m 2 (l + to 2 ). For 7 = 0, 1, 2 the 
Jeans equations can be solved explicitly, even in the presence of a central black hole, with 
functions no more complicated than standard hypergeometric functions. 

Appendix B: Relation with expansions in spherical harmonics 

Equation Q can be recovered also from a study based on spherical harmonics. 
Let us assume that the density distribution is given in spherical coordinates by 
p = — po rp'(r) cos 2 In a standard expansion in spherical harmonics p(x) = 
Sz^o Di m (r)Yi m (-d, if) (e.g., Jackson 1999); for our case we thus have 



AmW = -porp'{r)— ^-5 n 



610 + 7t Sl2 



(B.l) 



Some algebra then leads to prove that the associated potential is indeed the one recorded 
in Eq. J7J). 

In the seminal paper describing the numerical set-up of the "Schwarzschild method" 
for the construction of triaxial models of elliptical galaxies (Schwarzschild 1979), the 
(dimensionless) density distribution was chosen of the form 

Ps = F(r) G{r) 2z2 -£ 2 ~ y2 + 3H(r)*-=^, (B.2) 

where F(r) = (1 + r 2 )^' i / 2 , the two non-spherical terms have factors of the form of I = 2 
spherical harmonics, and H and G are functions prescribed implicitly from separate 
considerations and then calculated numerically. Given the above choice of p§, the related 
potential can be written as 

0s - U(r) V(r) 2z2 ~^ 2 ~ y2 + 3W(r)?-^-, (B.3) 

and then Schwarzschild obtained V = V[G] and W — W[H] by solving the (dimension- 
less) Poisson equation V 2 (/>s = ps numerically. Here we note that V[G] and W[fT] are also 
available in closed integral representation for a generic F(r) (starting from Eqs. [8]- [9] 
with the condition a + (3 + 7 = 0) . 

Several authors constructed triaxial density-potential pairs starting from Eqs. (B2)- 
(B3), through a potential-priority variation of the Schwarzschild method (e.g., de Zeeuw 
& Merritt 1983, Hernquist & Quinn 1989, de Zeeuw & Carollo 1996). We wish to em- 
phasize that the approach based on homeoidal expansions developed in this paper to 
construct analytically tractable density-potential pairs leads to new solutions. For exam- 
ple, models constructed from an expansion of triaxial 7-models are different from those 
studied by de Zeeuw & Carollo (1996). 
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Appendix C: Solutions of the Jeans equations for low-flattening, scale-free 
axisymmetric models 



Here we record the expressions for the intrinsic and projected kinematical profiles asso- 
ciated with the density profile in Eq. I|26|) . derived from the Jeans equations under the 
assumption that the underlying DF depends on E and J z . We adopt the same normal- 
ization procedure mentioned at the beginning of Sect. 4. 

We first note that for 7 > 1 the edge-on projected density ("surface brightness") is 



E = 



7^r( 7 /2 - 1/2) 

r(7/2)l7-i 



1 - 



(7 - l) r)z 



-72 



(C.l) 



vV 



and the mass contained within the projected radius t 
ity") is 

M _ tt 3 / 2 (2 + 77-77^(7/2- 1/2) - 3 _ 7 
(3- 7 )r( 7 /2) 

Integration of Eq. (|14|) . truncated to first order in 77, gives 



z 2 ( "integrated luminos- 



(C.2) 



go 



(5- 7 )(3- 7 )f 2 (-r- 1 ) 



Mbh 



5 — 7 — 477 (4 — 7)t7z 2 
2(7-1) 



2777 777Z 2 



(3 + 7)f 1 +T 
Similarly, from Eq. H15[l. 

— 2 2t7^ 2 - 
Q( v l - & ) - 



1 + 7 



1 



7MBH 



(C.3) 



(C.4) 



Therefore, for spherical systems (77 
order in 77 

ky/2rjR I 



(5-7)(3-7) (3 + 7)f 3 -T_ 

0), the quantity above vanishes. Finally, to leading 



Qu v = \ QX gu 



1 



7MBH 



f3 7 /2 



(C.5) 



(5-7)(3-7) (3 + 7 )f 3 -+ 
The relevant projected kinematical quantities are obtained from the formulae described 
in Sect. 3. Accordingly, 

2 V k 2 ^y 2 \ r( 7 -l/2) 7MBH r( 7 /2+l) 



£2 7 -l 

and, for 7 > 3/2, 

v^r( 7 



(5 - 7) (3 - 7 )r( 7 ) 2r( 7 /2 + 5/2) 



(C.6) 



1/2) 



(5-7)(3-7)r(7)^" 3 
/JBH^r( 7 /2 +1) /3 + 7 - 277; 



5 — 7 — 477 (4 — 7)77^2 
27-3 



£ 2 



777^2 



1 - k 2 
1 k EV 2 

k 2 p 



(C.7) 



2r(7/2 + 5/2)^7 V 7 e 2 
Unfortunately, when /ibh / no explicit expression is available for the quantities Eu p 
and Eitp. This latter quantity is obtained by expanding to lowest order in 77 the function 
(Eu p )2/E. In two interesting cases, namely when /jbh = and when /ibh 7^ and £ — > 0, 
asymptotic expressions can be obtained from the projection integral by substituting 



and then by changing the integration variable to r/£. One then finds 



/2nrjky 



r(3 7 /4-l/2) 



J t W-i v /(5- 7 )(3-7)r(3 7 /4)' 
^ V^iiTfcy VtT(7/2 + 1/4) 



7J7+1/2 V3T^r( 7 /2 + 3/4) ' 



(/Jbh = 0); 



(mbh ^ 0, 



(C.8) 
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and 




' 2^ V k 2 y 2 r 2 (3 7 /4-l/2)r( 7 /2) 

fa-i (5 - 7) (3 - 7 )r 2 (3 7 /4)r( 7 /2 - 1/2) ' 
2^ WBH k 2 y 2 7 r 2 ( 7 /2 + l/4)r( 7 /2) 

£ 2+7 (3 + 7 )r 2 ( 7 /2 + 3/4)r( 7 /2-l/2)' 



Obh = 0); 



(mbh ^0, *-0). 



(C.9) 



Near the origin, the black hole contribution to the various intrinsic and projected fields 
dominates for any value of 7 . 
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